function [Ahat, Resid, Sigmahat, Ahat_determ, X_determ, T_used, df_lost] = ...
               zFC_VAR_Estimation_exog(Y,p,constantD,trendD,trendBreaks,seasonD,moreBreaks,exogD,Y_exog)

% -----------------------------------------------------%
% DESCRIPTION OF THE FUNCTION
% this function estimates the reduced-form VAR model for the variables
% included in the matrix Y. 

% More precisely, the variables included in Y enter as rows
% variables. p lags are used for the endogenous variables. A constant is
% added if the dummy constantD equals 1. A trend is added if the dummy
% trendD equals 1. The trend is linear, unless quadraticD or cubicD
% indicate that the trend should be respectively quadratic or cubic. If
% timeBreaks is bigger than zero (and less than 2), then just as many break
% points are included in correspondence of the period(s) indicated in
% whenBreaks. An extra dummy variable is included in whenExtra if
% extraBreak equals 1. Seasonal dummies are included if seadonD equals 1.
% -----------------------------------------------------%
            
% check on stata 

% y_t = Ahat*[y_t-1; y_t-2; .... ] + Ahat_determ*X_determ_t + r_t

k      = size(Y,1);
T      = size(Y,2);
T_used = T-p;
assert(T > k) % check that that data at time t enter as column vector of Y

% Operative variables
y = Y(:,p+1:T)'; 
X = Y(:,p:T-1)';
for g = 1:p-1
   X = [X, Y(:,p-g:T-1-g)'];  
end

% constant
if constantD == 1
    X = [X, ones(T-p,1)];
end

% exogenous variables
if exogD == 1
   X = [X, Y_exog(:,1+p:T)']; 
end

% time trend
if trendD >0 
    step = [1:1:size(X,1)]';
        if trendD == 2;
            step = step.^2;
        end
        if trendD == 3;
            step = step.^3;
        end
        step = log(step); % without log the non linear trends give a X that is hard to invert
        X = [X, step];  
    
        % add break points only if dataset includes the corresponding time
        % breaks. Otherwise, adjust break dummy so that break is not
        % included
        
        % rewrite so that it is no proble if breaks are not given or if
        % they are just a NaN


        if length(trendBreaks) > 0
             time_step = [1:1:T]'; % observation in the original Y where we add the break
             assert(length(time_step) == size(Y,2))
                 for ttt = 1:length(trendBreaks)
                     step2 = step;
                     step2(1:trendBreaks(ttt)-p) = zeros(trendBreaks(ttt)-p,1);
                     X = [X, step2];  
                 end
                 % check.........
        end

end

% seasonal dummies
if seasonD == 1
    step = eye(4);
    if constantD == 0
    else
        step = step(:,2:4);        % to avoid dummy variable trap
    end
    step = kron(ones(floor(T/4)+2,1),step);  
    step = step(p+1:T,:);  
    X = [X, step];
    step_season = step;
end

% an additional break points
if moreBreaks > 0
       time_step = [1:1:T]'; % observation in the original Y where we add the break
       assert(length(time_step) == size(Y,2))
       for ttt = 1:length(moreBreaks)
           step = zeros(size(X,1),1);
           step(moreBreaks(ttt) - p) = 1;
           X = [X, step];  
       end
else
end

df_lost = size(X,2);
T_used = T-p;

% save X variable, so that it is easy for the bootstrap. and summarize
% deterministib bits into a constant?


% ESTIMATE THE MODEL
% Ahat_step = (X'*X)^(-1)*X'*y;     
 Ahat_step = (X'*X)\(X'*y);     

Resid = y-X*Ahat_step;
RSS = Resid'*Resid;                       
Sigmahat = RSS/(T_used-df_lost);    % variance covariance of the reduced-form shocks

% EXPORT ESTIMATES
Ahat = Ahat_step(1:k*p,:)';
Resid = Resid';

if size(X,2) > k*p
    X_determ = X(:,k*p+1:end);
    Ahat_determ = Ahat_step(k*p+1:end,:)';
else
    X_determ = 0;
    Ahat_determ = zeros(k,1);
end


end